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ABSTRACT 


When  a  plane  harmonic  wave  has  oblique  incidence  on  a  stack  of 
homogeneous  dielectric  slabs,  the  field  in  each  layer  consists  of  an 
outgoing  and  a  reflected  wave.  The  complex  amplitudes  of  these  plane 
waves  in  layer  n  are  denoted  by  An  and  Bn,  and  simple  recursive  re¬ 
lations  are  derived  for  An+l  and  Bn+i  as  functions  of  An  and  Bn.  The 
solution  begins  by  setting  AQ  =  1  and  BQ  =  0,  and  it  is  completed  by 
calculating  the  transmission  and  reflection  coefficients  of  the  plane 
multilayer  as  simple  functions  of  Aj^+j  and  Bj^+j  where  N  represents 
the  total  number  of  layers. 

This  technique  is  also  applicable  when  a  plane  wave  is  incident 
on  a  multilayer  dielectric  cylinder  or  sphere.  The  required  equations 
are  derived  for  the  cylindrical  case  with  normal  incidence.  The  rigorous 
solutions  and  digital-computer  programs  are  included  for  both  principal 
polarizations  for  plane  and  cylindrical  multilayers. 
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EFFICIENT  RECURSIVE  SOLUTIONS  FOR 
PLANE  AND  CYLINDRICAL  MULTILAYERS 


I.  INTRODUCTION 

The  transmission  and  reflection  coefficients  of  a  plane  dielectric 
multilayer  are  pertinent  in  the  design  of  radomes  and  in  other  appli¬ 
cations.  One  type  of  corner  reflector  has  one  of  its  three  conducting 
surfaces  coated  with  several  dielectric  layers  which  arc  designed  to 
produce  a  right -circularly  polarized  reflection  when  the  incident  wave 
is  right -circularly  polarized.  Multilayer  dielectric  spheres  and  cy¬ 
linders  form  an  important  class  of  antenna  scanning  and  echo  enhance¬ 
ment  devices. 

A  matrix -multiplication  solution  for  the  plane  multilayer  is 
described  by  Collin[l].  This  matrix  solution  represents  a  distinct 
advance  over  the  older  methods.  A  recursive  formulation  developed 
here  is  even  more  straightforward  in  its  theory  and  more  efficient  for 
numerical  calculations.  The  technique  is  applicable  also  to  radially 
layered  cylinders  and  spheres,  waveguides  containing  two  or  more  media, 
and  surface  waves  on  plane  multilayers. 

This  report  presents  the  derivation  of  the  recursion  formulas  for 
plane  and  cylindrical  multilayers.  The  two  principal  polarizations  are 
considered,  and  digital -computer  programs  are  included.  This  solution 
can  be  applied  directly  to  the  analyses  of  inhomogeneous  slabs  and  cy¬ 
linders  by  using  piecewise  uniform  approximations  for  the  permeability 
and  permittivity  functions. 


II.  THE  THEORY  OF  THE 

PLANE  MULTILAYER 

Suppose  a  harmonic  plane  wave  in  free  space  has  oblique  incidence 
on  a  plane  multilayer  consisting  of  N  homogeneous  isotropic  slabs,  as 
indicated  in  Fig.  1.  Lei  d^,  ja.n  and  cn  represent  the  thickness,  permea¬ 
bility,  and  permittivity  of  slab  n.  The  slabs  are  considered  to  have 
infinite  width  and  height  and  parallel  surfaces,  with  unbounded  free 
cpace  on  both  sides  of  the  multilayer.  The  incident  plane  wave  impinging 
on  the  left-hand  surface  of  the  multilayer  is  given,  in  the  TE  case  (i.  e. 
perpendicular  polarization)  by 
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Fig.  1.  A  plane  multilayer,  illustrating  the  outgoing  and 
reflected  waves  in  each  layer. 


(1) 


ejkDy  sin6  ejkGz  cos9 


where  0  is  the  angle  of  incidence*  kD  =  2tt/X,  and  X  is  the  free- space 
wavelength.  The  reflected  plane  wave  is  given  by 


(2)  EJ  =  R  Eq  eik°y  sin0  e-*»koz  CO80  , 


where  R  is  the  reflection  coefficient  of  the  multilayer.  The  transmitted 
plane  wave  on  the  right-hand  side  of  the  multilayer  is  represented  by 


(3) 


E‘  =  T  Eq  ejkoV  sin0  eJkoz  cos°  , 


where  T  is  the  transmission  coefficient  of  the  multilayer.  The  field  in 
each  layer  can  be  regarded  as  an  infinite  series  of  plane  waves  bouncing 
back  and  forth,  but  it  is  more  convenient  (and  equally  valid)  to  consider 
it  to  be  the  sum  of  only  two  plane  waves,  one  traveling  outward  and  one 
reflected.  v 
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In  layer  n,  for  example,  the  field  is  represented  by 


(4)  En  =  (An  e^2  +  Bn  e’Yn  =■  . 

Similarly,  in  layer  n+1  the  electric  field  intensity  is  given  by 

(5)  En+1  =  eVn+lJ  t  Bn+1  e^n+1  ejk0y  . 

The  boundary  between  layers  n  and  n+1  is  located  at 

(6)  zn  =  <*i  +  ^2  +  ds  +  +  +  dn 

By  enforcing  the  boundary  conditions  on  Ex  and  Hy  at  z  =  zn,  we  can  show 
that 

(7)  An+i  =  Pn  An  4  Q,,  Bn 
and 

(8)  Bn+1  =  Bn  ^xi  +  Sn  Bn  , 
where 

(9)  Pn  =  °.  5  (1  +  yn+l  Yn/^n+l)  e(Yn-Yn+]  >  zn  , 

(10)  On  =  0.5  (1  -  4n+iYn^nYn+1>  e'^n*  Yn+1>  z„  # 

(11)  Rn  =  °.  ^  (i  -  &n+lYn^nYn+i)  elYn  +  Yn+l)zn  # 

and 

(12)  Sn  =  0.5  (1  +  |ln+1Yn/^nYn+1)  e‘(Yn  -  Vntl^n  % 
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The  propagation  constant  Yn  for  layer  n  will  be  complex  if  the 
medium  is  dissipative.  Both  the  real  and  imaginary  parts  of  Yn  will  be 
positive.  If  layer  n  is  a  lossless  medium,  Yn  will  be  purely  imaginary. 
The  wave  equation  is  employed  to  obtain  the  following  equation  for  yn: 


(13)  Yn  =  J  An<n  "  ko  sin2  0  * 


The  reflection  and  transmission  coefficients  of  the  multilayer 
(R  and  T)  are  often  of  particular  interest.  We  can  calculate  these 
quantities  in  a  systematic  manner  by  setting 


(14)  Aq=1 
and 

(15)  BQ  -  0 


and  then  using  the  recursion  equations  Eqs .  (7)  and  (8)  to  calculate  Aj, 
Bi,  A2,  B*,  ...  Aivf+1,  and  Bn+1  *n  that  order.  This  process  is  usually 
carried  out  on  a  digital  computer. 

From  Eqs.  (1)  through  (4), 


(16) 

Eo  =  AN+1  * 

(17) 

R  _  bN+1  _ 

an+i 

and 

(18) 

T  =  1 

AN+1 

In  the  TE  case  the  constants  An  and  Bn  represent  the  electric 
field  intensities  of  the  outgoing  and  reflected  waves  in  each  layer.  In 
the  TM  case  (parallel  polarization)  the  solution  proceeds  in  the  same 
manner.  In  fact,  the  equations  given  above  apply  in  both  cases  but  the 
An  and  represent  the  magnetic  field  intensities  in  the  TM  case  and 
P-n+l  and  pn  must  be  replaced  with  «n+l  and  «n  in  Eqs.  (9)  through  (12). 

If  a  perfectly  conducting  sheet  is  placed  on  the  right-hand  surface 
of  the  multilayer  (i.  e. ,  on  the  xy  plane),  the  solution  is  again  given  by 
the  equations  above  with  the  exception  that  we  do  not  calculate  the  trans¬ 
mission  coefficient  T  in  this  case,  and  Eqs.  (14)  and  (15)  are  replaced 
with 


(19)  Aj  =  1  and  Bj  =  -1  in  the  TE  case 
and 

(20)  Aj  =  1  and  Bj  =  1  in  the  TM  case. 


Equations  (19)  and  (20)  are  obtained  by  forcing  the  tangential  electric 
field  intensity  to  vanish  at  the  perfectly  conducting  plane. 

The  insertion  phase  delay,  denoted  by  the  symbol  IPD,  is  equal  to 
the  phase  angle  of  the  complex  transmission  coefficient  T  but  is  of  opposite 
sign: 


(21)  IPD  =  -  Phase  (T) 

In  the  previous  equations,  the  reflection  coefficient  R  is  defined  as 
the  ratio  of  the  reflected  wave  amplitude  to  the  incident  wave  amplitude 
at  the  coordinate  origin;  that  is, 

E*  (0,0,0) 

(22)  R  -  — - -  for  the  TE  case 

0,0) 

and 

h!(o.  o,o) 

(23)  R  =  — - -  for  the  TM  case. 

*4(0,  0,0) 
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Let  R  represent  the  ratio  of  the  reflected  and  incident  wave  amplitudes 
at  the  point  of  incidence  where  (x.  y,  z)  =  (0,  0,  d),  with  d  representing  the 
total  thickness  of  the  multilayer.  Then 


(24)  R'  =  R  e -2jkod  cos  6. 


Equation  (24)  can  be  verified  by  means  of  Eqs.  (1)  and  (2)  and  the  de¬ 
finition  of  R'.  The  two  reflection  coefficients,  R  and  R',  differ  only  in 
phase. 

A  digital- computer  program  based  on  these  equations  is  given  in 
Appendix  I,  and  some  numerical  results  are  included  in  Section  IV. 


III.  THE  THEORY  OF  THE  CYLINDRICAL 
MULTILAYER 

Consider  a  plane  harmonic  wave  in  free  space  to  have  normal 
incidence  on  a  dielectric  cylinder  of  infinite  length,  as  suggested  in  Fig. 
2.  The  cylinder  axis  is  taken  to  be  the  z  axis  in  a  rectangular  coordinate 
system,  and  the  x  axis  is  the  axis  of  propagation  of  the  incident  plane 
re.  Let  the  cylinder  consist  of  M  lossless  homogeneous  layers,  each 


PROPAGATION 
AXIS  * 


x 


Fig.  2.  Circular  dielectric  cylinder  with  several 
homogeneous  layers. 
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layer  being  a  circular  cylindrical  shell.  The  permeability  and  permit¬ 
tivity  of  layer  m  are  denoted  by  pm  and  «m,  and  the  phase  constant  is 
given  by 


(25)  km  =  uJVm  *m  • 


In  the  TM  case  the  electric  field  intensity  has  only  a  z  component,  given 
in  layer  m  by 


(26)  Em  =  ^  [Amn  Jn(kmp)  +  Bmn  Nn(kmp)J  cos  n<p  , 

n=0 


where  (p,$,  z)  are  the  cylindrical  coordinates  and  JR  and  Nn  represent 
the  Bessel  and  Neumann  functions.  The  time  dependence,  ej^t,  is  under¬ 
stood.  In  the  TE  case  the  solution  is  obtained  from  the  equations  given 
here  by  interchanging  p  and  c  and  E  and  H,  where  H  represents  the  mag¬ 
netic  field  intensity. 

The  coefficients  Amn  and  Bmn  must  be  determined  by  applying  the 
boundary  conditions  on  Ez  and  at  each  interface.  From  Maxwell's 
equations,  the  component  of  the  magnetic  field  intensity  in  layer  m  is 
given  by 


(27)  Hm  =  (km/jwpm)  Y  [Amn  JR(kmp)  +  Bmn  N^(kmp)J  cos  n«|>  , 

n=0 


where  JR  and  denote  the  derivatives  of  the  Bessel  and  Neumann  functions 
with  respect  to  the  argument. 

Let  Rm  denote  the  outer  radius  of  layer  m.  From  the  boundary 
conditions  at  this  interface  between  layers  m  and  m+1,  it  is  found  that  the 
wave  expansion  coefficients  for  the  two  regions  have  the  following  linear 
relations: 


(28) 

and 

(29) 


Am+l,n  =  UmnAmn  +  WmnBmn 


®m+l,n  =  ^mn^mn  +  ^mn®mn  * 
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where 


*  ^m^rn+lJn^m^m^n^m+l^-rn^  * 

(31)  vmn  "  (TrRm^|im)  [  ^m+ l^mJn^mRm)Jn^kmT  1  Rm) 

^m^m + 1  ^n^m+ 1  Rm  ^J  » 

<32>  Wmn  -  ^+lkmN^kmRm>Nn<km+lRm> 

+  ^mkm+lNn(kmRm)N;(kmilRm)J  . 

and 

(33)  ^mn  “  [  ^m+l^m^n^m^m^n^m+l^m^ 

-  u  k  ,,N  (k  R  )j'(k  ±,k  )1  . 

m+1  n  m  m  n  m+1  m  J 


If  the  coefficients  in  the  first  layer  (Ajn  and  Bin)  were  known,  the 
coefficients  in  the  remaining  layers  could  be  calculated  by  using  Eqs.  (28) 
and  (29)  recursively.  To  permit  a  procedure  of  this  type,  let  us  define  a 
set  of  normalized  coefficients  A^nn  and  B^^  related  to  Amn  and  Bmn  by 
the  proportionality  constant  Kn  as  follows: 


(34)  Amn  =  Kn  A^^ 

and 

(35)  Bmn  =  Kn  Bmn  • 


With  no  loss  of  generality,  let 


(36)  Ajn  =  1  . 
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If  the  center  layer  is  a  dielectric  medium,  the  field  must  be  finite  at 
the  origin  and 


(37)  =  0  . 


The  normalized  coefficients  also  obey  the  same  recursion  formulas  and 
we  can  now  calculate  A]n,  B^r,  ...  A^n  ,  B^ln.  AM+l,n*  and  BM+l,n 
in  that  order. 

The  field  in  the  exterior  free -space  region  is  given  by 

o* 

(38)  Em+1  =  Y  [(-j)n  en  Jn(koP)  +  CnHn(!,(koP)]  cos  no  . 

n=G 

where  Hn^(k0p)  represents  the  Hankel  function  and  kQ  is  the  phase  con¬ 
stant  of  free  space.  The  first  series  in  Eq.  (38)  represents  the  incident 
plane-wave  field,  and  the  second  series  is  the  scattered  field  which  con¬ 
tains  outward-traveling  waves  only.  The  function  en  is  unity  if  n  =  0,  and 
en  =  2  if  n  is  greater  than  zero. 

The  field  in  the  exterior  free-space  region  is  also  given  by  Eq.  (26) 
with  m  =  M+l.  Comparison  of  the  two  representations  reveals  that 


(39)  AM+l,n  -  J  BM+l,n  =  (-j^  en  * 


The  constant  of  proportionality,  Kn,  is  found  from  Eqs.  (34),  (35),  and 
(39)  to  be 


(40) 


(-j)n 


Kn=  -r 


'n 


AM+1,  n"J  BM+l,n 


The  scattering  coefficients  for  the  external  region  are  given  by 


(41) 


Cn  =  j  B 


H)n  en  Bi 


M+l,n 


M+l.n 


B' 


+  j  Av 


M+l,  n  J  M+l,  n 
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This  completes  the  solution.  Equation  (26)  can  be  employed  to  calculate 
the  field  at  any  point  in  the  dielectric  cylinder,  and  Eqs.  (38)  and  (41) 
are  used  to  calculate  the  field  at  any  point  outside  the  cylinder. 

If  the  dielectric  media  are  lossless,  it  may  be  noted  that  the  con¬ 
stants  A^n,  B^n,  Umn,  Vmn,  Wmn.  and  Xmn  are  real.  The  external 
scattering  coefficients,  Cn,  are  complex. 

In  the  case  of  a  perfectly  conducting  circular  cylinder  with  one 
or  more  dielectric  layers,  let  "a"  be  the  radius  of  the  conducting 
cylinder  and  Rj  the  outer  radius  of  the  first  dielectric  layer.  The  above 
equations  again  give  the  solution  if  Eq.  (37)  is  replaced  with 

Jn(kia) 

(42)  B{n  =  -  -  for  the  TM  case, 

Nn(ki  a) 

and 

,  J' (ki  a) 

(43)  Bin  =  ~  — - _ for  the  TE  case. 

Nn(kia) 

Equations  (42)  and  (43)  are  obtained  by  setting  Ez  =  0  or  Eq>  =  0  at  the 
conducting  surface. 

At  any  point  in  space  outside  the  dielectric  cylinder,  the  scattered 
field  is  given  in  the  TM  case  by 

oa 

(44)  Es  =  )  (~j)n  en  DnHn(k0p)  cos  n4>  , 

n=0 


where  the  superscript  (2)  is  understood  on  the  Hankel  functions  H  (kQp)  and 


(45) 


~BM+1,  n _ 

BM+ 1,  n  +  j  AM+l,n 


In  calculating  the  scattered  field  at  a  great  distance  from  the  cylinder,  we 
use  the  asymptotic  form  for  the  Hankel  functions  of  large  argument  to  show 
that 


(46)  Es  =  V2j/ TTk0p  e’^kot>  ^  enDn  cos  n<?  • 

n=0 
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When  a  plane  wave  is  incident  on  a  cylindrical  structure  of  infinite 
length,  the  distant  scattering  pattern  is  conveniently  described  by  the 
echo  width  which  is  defined  as  follows: 


(47)  W  =  limit  2frp  ,ES/Ei'2  . 

p  —  « 


In  Eq.  (38)  the  incident  electric  field  intensity  El  is  taken  to  have  unit 
magnitude.  From  Eqs.  (46)  and  (47)  the  bistatic  echo  width  of  the  multi¬ 
layer  dielectric  cylinder  is  given  by 


(48) 


W  =  (2\/ir) 


y 


n=0 


enDn  cos 


where  X.  is  the  wavelength  in  free  space. 

A  digital-computer  program  based  on  these  equations  is  included 
in  Appendix  II. 


IV.  NUMERICAL  RESULTS 

Figure  3  shows  the  transmission  coefficient  and  insertion  phase 
delay  versus  frequency  for  a  plane  multilayer*  calculated  with  the  aid 
of  the  equations  in  Section  II.  Numerical  calculations  obtained  with  these 
equations  have  been  found  to  agree  with  those  published  by  other  investi¬ 
gators. 

Figure  4  ~hows  the  distant  scattering  pattern  of  a  cylindrical 
multilayer,  calculated  with  the  equations  given  in  Section  III.  These 
equations  are  also  found  to  yield  results  which  agree  with  previously 
published  data. 


V.  CONCLUSIONS 

An  efficient  recursive  solution  is  developed  for  the  transmission 
and  reflection  coefficients  of  a  plane  multilayer  and  for  the  scattering 
pattern  of  a  cylindrical  multilayer.  The  two  principal  polarizations  are 
considered,  and  digital-computer  programs  are  included.  The  technique 
is  also  applicable  to  the  multilayered  sphere. 
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The  appropriate  equations  are  also  given  for  the  reflection  co¬ 
efficients  of  a  perfectly  conducting  plane  which  is  coated  with  a  stack 
of  homogeneous  dielectric  sheets,  and  for  the  scattering  pattern  of  a 
perfectly  conducting  circular  cylinder  coated  with  several  homogeneous 
dielectric  layers. 
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Calculated  power  transmission  coefficient  and  insertion 
phase  delay  versus  frequency  for  plane  multilayer 
having  three  layers. 


IPO  (DEGREES) 


ECHO  WIDTH  /WAVELENGTH 


<fi  (DEGREES) 


Fig.  4.  Distant  scattering  patterns  of  a  five-layer  dielectric 
cylinder,  calculated  with  the  recursion  technique. 


14 


APPENDIX  I 

COMPUTER  PROGRAM  FOR  THE  PLANE  MULTILAYER 


A  computer  program  for  the  reflection  and  transmission  coefficients 
of  a  plane  multilayer  is  shown  in  Fig.  5.  This  program,  written  in  the 
computer  language  known  as  Scatran,  is  based  on  the  equations  given  in 
Section  II.  The  symbols  used  for  the  input  data  are  defined  as  follows: 


KKK 


N 

D(I) 

E(I) 

TD{I) 

U(I) 

TP(I) 

THETA 

FGC 


=  Number  of  cases  to  be  calculated  with  different 
frequencies  or  angles  of  incidence 
=  Number  of  layers 
=  Thickness  of  layer  i  in  inches 
=  Dielectric  constant  of  layer  i  relative  to  that  of 
free  space 

=  Electric  loss  tangent  of  layer  i 

=  Permeability  of  layer  i  relative  to  that  of  free  space 
=  Magnetic  loss  tangent  of  layer  i 
=  Angle  of  incidence  in  degrees 
=  Frequency  in  gigacycles 


Some  of  the  other  symbols  used  in  the  program  are  defined  below. 


TH 

WAVE 

SS 

CC 

UC(I) 

EC(I) 

G(  I) 

Z(I) 


=  Angle  of  incidence  in  radians 
=  Wavelength  in  free  space,  inches 
=  sin2  G 
=  cos  9 

=  Complex  relative  permeability  of  layer  i 
=  Complex  relative  permittivity  of  layer  i 
=  y/kQ  for  layer  i,  as  given  by  Eq.  (13) 

=  kQz{  where  z^  is  defined  by  Eq.  (6) 


The  field  intensities,  An  and  Bn,  of  the  outgoing  and  reflected  waves  in 
layer  n  are  denoted  by 


AE  and  BE  for  the  TE  case  without  a  conducting  plane, 
AM  and  BM  for  the  TM  case  without  a  conducting  plane, 
AEC  and  BEC  for  the  TE  case  with  conducting  plane,  and 
AMC  and  BMC  for  the  TM  case  with  conducting  plane. 


Between  statements  S30  and  S50,  the  recursive  calculations  are  carried 
out  in  accordance  with  Eqs.  (7)  through  (12).  Immediately  following 
statement  S50,  the  transmission  and  reflection  coefficients  are  calculated 
by  means  of  Eqs.  (17),  (18),  and  (34).  These  coefficients  are  denoted  by 
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TE 

TM 

RE 

Rivl 

REC 

RMC 


=  Transmission  coefficient  for  the  TE  case  without 
conducting  plane, 

=  Transmission  coefficient  for  the  TM  case  without 
conducting  plane, 

=  Re  flee  tion  coefficient  R1  for  the  TE  case  without 
conducting  plane, 

=  Reflection  coefficient  R1  for  the  TM  case  without 
conducting  plane, 

=  Reflection  coefficient  R'  for  the  TE  case  with  a 
conducting  plane,  and 

=  Reflection  coefficient  R1  for  the  TM  case  with  a 
conducting  plane. 


The  symbols  used  for  the  output  data  are  defined  by 


REP,  PRE,  TEP,  FIPDE 

RMP,  PRM,  TMP,  FIPDM 

RECP,  PREC 

RMCP,  PRMC 


=  Power  reflection  coefficient, 

reflection  phase,  power  transmission 
coefficient,  and  insertion  phase  delay 
for  TE  case  without  conducting  plane; 
=  Power  reflection  coefficient,  reflec¬ 
tion  phase,  power  transmission  co¬ 
efficient,  and  insertion  phase  delay 
for  TM  case  without  conducting  plane; 
=  Power  reflection  coefficient  and 
reflection  phase  for  TE  case  with 
conducting  plane;  and 
=  Power  reflection  coefficient  and 
reflection  phase  for  TM  case  with 
conducting  plane. 


A  typical  set  of  input  data  in  the  proper  format  is  shown  at  the  end  of  the 
computer  program. 

This  program  will  handle  a  maximum  of  100  layers,  but  this  number 
can  be  increased  to  a  much  larger  number  simply  by  modifying  the  dimen¬ 
sion  statement  near  the  beginning  of  the  program.  Many  of  the  statements 
can  be  deleted  when  the  multilayer  with  a  conducting  plane  is  the  only  case 
of  interest,  or  when  this  case  is  of  no  interest.  Furthermore,  the  program 
can  be  simplified  if  all  of  the  layers  have  the  same  permeability  as  free 
space.  A  simple  modification  in  the  program  can  provide  for  increments 
to  be  taken  in  the  frequency  or  the  angle  of  incidence. 
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*■*# 

#** 


FIRST 


S<5 


520 


530 


RON 

SCATRAN 

COMPLFX (C5QRTL. .CEXPL. .  C  I  « t  tC  ,  FO  ,  G  .  A^  ♦  A M  ,  qr  ,  rm  ,  AFC  %  Rrr  ,  07  % 
AMC.BMC.F.YI  ,  VI  I  «  EXP  1  «  EXP  !  I  ,  AEP.  AMP»ryp,pr,PM,AFCD«PfrC« 
AVCP.RVO.FP. FV.TS.TM)- 

DIMFNSI ON(D(  1  ''O  >  ,(T(  J  )  ,  TD  (1  OO  )  ,1)  (  1  )  ,  TP(  I  no  )  .7  (  ir>n  )  , 

UCC 1 ) ,EC( 1 n'M  ,G(  1  00)  >- 
READ  1NPUT.S. (<K<«N)— 

R^  AD  INPUT. 7,  (  (D(  I  )  «  I  =  1  .  1  .  !  .|_F.M)  )- 
RFAD  INPUT.  7.  (  (  R  (  1  >.1*1  .1  «1*LP«N)  )  — 

RFAD  INPUT. 7.  (  (TO(  !  ). 1  =  1*1  «  !  .Lr • N )  )  - 
READ  INPUT.7,HlMh«!s|.l.!.LF«N)l- 
RFAD  INPUT. 7.  (  ( TP (  I  ) « 1 = 1 « 1  .  I  .LF.N  )  >- 
NN=  N+ 1 - 

no  through (ss i « i = i « i . i .le.n- 

F  I  =  I  — 

WRITE  OUTPUT. 2* (FI«D(I ) «F( 1 ) ,TD( I )  .il(  1 )  .  TP( I ))- 

TP  I  =6  •  283 1  85-*- 

DR=B7«  295779- 

CI«.I«1.- 

c'CCNNIM.- 

UC (NN)=1 •- 

Z ( 0 ) = • 0- 

WRITF  OUTPUT. 2- 

DO  THR0UGH(  51  00  )  .<=  1  .1  «K*LE  .«K- 

read  input«7. (theta, fgc i- 

TH  =  «  0 1 74 5 329* T HE T  A  — 

WAVE  =  1 1 .800! 4/FGC- 

WRITE  OUTPUT,  2.  (  T  HP  T  A  »  F  GO  ,  A  VF  )- 

WRITE  OUTPUT, 2- 

C5=5IN. ( TH)*5 IN. ( th )- 

CC=COS. (TH)- 

OO  THROUGH!  O?''  )  ,  I  =  1  .  I  ,  1  .Lp  .N- 

uc  ( i  ism  i  )*(i.-ci#TPcn*- 
pr ( i ) =f ( i ) * ( i .-r i*td ( 1 1 ) - 

G(  I  )  =  C  I  *C 50RTL  •  (Uf{  1  )*rn  I  )-cc  I- 

gz=c i *cr- 

G  (NN  )  =C  1  *C.C- 

no  through (57^)»i*i»i»i.le ,n- 
11*1-1- 

Z ( I  ) =Z (  I  I  ) +D (  I  )*TPl/WAVr- 
AF= . 5* (  1 .+ur (  1  )*G7/G ( 1  )  >- 
PFs . R# (  1 .-urc 1  )#G Z/Gf  I’l- 
AV=.R*{ l  .+er  (  1  )*0Z/C  (  1  )  )- 

pv=.r*  ( i  .-er  ( i  )*g 7/r- ;  r  •;  )- 
AFC  = 1 •  — 

BFC  =- \ •- 
AMC= 1 • - 
RVC  =  1 •- 

DO  THROUGH ( SR"'  )  »  I  I  =2 .  1  ♦  I  I *LE .NN- 
I  =  I  I -1  - 

F  =UC (II  )  *G ( 1  ) / ( UC (  I  ) *G (11))- 
Y I =  G (  I  )*Z(  1  )- 
Y  I  I  =G (  I  1  ) *7(1)- 


Fig.  5.  Digital- computer  program  for  reflection  and  transmission 
coefficients  of  plane  multilayer. 
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550 


5100 

*** 

1 


r p= i • +f_ 

F v= 1 *-r_ 

FxP  1  =CFX=>L  #>'v!  )  - 

fxpi  Ur.PxPt.  <  yi  I  i- 

A«=T>  =  AF- 
A  VP  =  A  V- 
Acr(-p=AF(-_ 

A  MCP=  A’JIO- 

A^=  (  .5/FXPJ  I  )KAfp«rxP|*c"'  +  '”  *PV/rxn|  )  , 
nr  =  .R*fxP  I  I  *  f  A^P*i- XP  !  *r  v  +  pt  *f;  >/r  yr  I  )  _ 

AE0  =  (  •F/f'xP  !  I  )  K  C  Arr  p»t-xp  (  *rr-4 •'<«-/-* fv/i  y p  ;  )  - 
PFCs»5*FXP  I  I  *  (  A-rr*'  xc  l  *'  v  +  pt-r  irp/r  yp  |  )_ 

F  =F r  (  1  I  )*G( I  >/(^C ( !  )*C (  I  I  )  )- 
FP= 1  • ♦  F  — 

Fv  = i . -F- 

A»s  (  .5 /FXP I  I  )  *  (  AVD*c7Xo  I  *rP  f(V/*crM/c-yr> ;  )_ 

P  vt=  *R*FxP  I  i  *  (  a  MP*c-yr>  i  #r  »/  +  ;-.«■» fo/p  yp  i  )  _ 

A VC  = ( .5/FxPI I  ) * < AVCP*F  XP I *r P+«vc*F v/(  yr i  )_ 
BMC  s  *5*FXP  I  I  *  (  AVCP*C  x°!  n"v+Mvc*r  P/”XP  !  )- 
TF=  1  •  /  A  F  -» 

TM= 1 ./AM- 

ARG=2**CC*Z(N)- 

FXP  =  C05.  ( ARG  >-• I  .  G|N.  (  A RG  ) - 

Rf=FXP*«^F/AF- 

qm=fxp*RM/Am_ 

RFC=FXP*BFr/AFC- 
Rw-C  =  CXP*R  MC/  A  Mf_ 

TFP= ( .AR5.TF ) .P.P- 
TVP= ( .ARS.TV ) .P.2- 
RFP  = ( . ARS.RE ) .P.2- 
RVP=(«AR5*RV).P*2- 
RFC.P  =  (  •  &=iS»Qr-C  )  .P.2- 
RMCP= ( • AHS.^vr  » .P.2- 

F  IP0F  =  -0R*FATAV?,  (  •  I  MAG.  T  "  «  .R-'  &L  «Tr  )  - 
f  I  Pnv  =  -0R*F  A  T  an?.  (  •  !  VA5.  rv«  •p-al.T’*)- 
PRFrDR*FA  T  AN2  •  (  •  I  V  AC.,Rt-  ,  .or  AL  «R~  »  - 
PRM  =  nR*F  A  T  AN 2 • < •  I  v  A  G  •  R  v  • .Rf.v_.Rv )_ 

PRFf  =0R*F  A  TAM?  •  (  •  iVA'.f?^,  A*_  ,Rr''!- 
PRMr  =  DR*FATA\i2.  (  •  !  VAG.Pm-  ,  . R " ,R\*'-  )  - 
WR  I  TF  0(  JTP<  JT  ,  2  «  <  Rc-P,  PRF  ,  T*  e*  f  |  .  _ 

WR  1  TF  Ot  ITPi  IT  «  2  «  (  *  PR  v  «  T'*o.K  ’  P-'V  j  . 

WRI  TE  O'  IT P', IT  «  2  «  C  RF  ~P  «  PPf  C.  >_ 

WRITE  OUTP'.  IT  «  2  «  (RvfO.PRvC  )- 
WRITE  OJTP'lT.2- 
FNO  PROGRAM ( F JRGI )- 
PATA 


•  I  5 

•  1  " 

•  1  R 

•  25 

*  OF 

c 

• 

u 

a  •  o 

2.0 

I 

^  p 

•  *  ^  2 

•  o-* 

.  *'  V 

1  •  ~ 

1  •  ~ 

1  •  - 

1  .  " 

1  .n 

•  w 

• 

•  ’ 

• 

• 

•  0 

io.- 

Fig.  5.  Digital-computer  program  for  reflection  and  transmission 
coefficients  of  plane  multilayer,  (cont.  ) 
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APPENDIX  II 

COMPUTER  PROGRAM  FOR  THE  CYLINDRICAL  MULTILAYER 


Figure  6  shows  a  computer  program  in  Scatran  for  calculating  the 
scattering  coefficients  and  the  scattering  pattern  of  a  layereo  dielectric 
cylinder,  based  on  the  equations  developed  in  Section  III.  The  input 
data  symbols  are  defined  as  follows: 

FN  =  Number  of  layers, 

DPH  =  Increment  in  scattering  angle  <p  in  degrees, 

RL(I)  =  Outer  radius  of  layer  i  in  free-space  wavelengths, 

E(I)  =  Permittivity  of  layer  i  relative  to  that  of  free  space,  and 

U(I)  =  Permeability  of  layer  i  relative  to  that  of  free  space. 

This  computer  program  is  applicable  only  to  lossless  dielectric  cylinders 
in  which  the  center  region  is  dielectric.  With  a  slight  modification  it  will 
give  the  solution  appropriate  for  a  perfectly-conducting  center  region. 

In  the  form  shown,  the  program  v.  Ill  handle  a  maximum  of  100  layers, 
but  this  can  be  extended  greatly  by  changing  the  dimension  statement. 

In  this  program  the  highest  order  of  the  Bessel  and  Neumann 
functions  required  is  determined  from  the  approximate  equation  L  =  4 
+  kj^Rjvj*  where  N  denotes  the  number  of  layers.  In  its  present  form, 
the  program  sets  an  upper  limit  of  100  for  the  maximum  order  of  the 
Bessel  and  Neumann  functions. 

The  wave  amplitudes  in  the  various  layers  are  denoted  by  A(I)  and 
B(I)  for  the  TM  case  and  AP(I)  and  BP(I)  for  the  TE  case,  where  I  re¬ 
presents  the  mode  number  n.  Between  statements  S25  and  S26,  the  program 
sets  A|n  =  1  and  Bjn  -  0  in  accordance  with  Eqs.  (36)  and  (37).  Next, 
the  arguments  kmRm  and  km+jRm  for  the  Bessel  functions  are  calcu¬ 
lated  and  are  denoted  by  W.  Between  statements  S26  and  S50,  the  Bessel 
and  Neumann  functions  are  obtained  by  calling  a  subroutine.  The  follow¬ 
ing  notation  is  used: 

FJA(I)  =  Jt(kmRm), 

FJB(I)  =  Ji(km+1Rm). 

FNA(I)  =  Ni(kmRm),  and 

FNB(I)  =  Ni(km+1Rm). 

In  its  present  form,  the  subroutine  for  the  Bessel  and  Neumann  functions 
is  limited  to  arguments  smaller  than  1000. 
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Between  statements  S50  and  S80,  the  derivatives  of  the  Bessel 
and  Neumann  functions  are  calculated.  These  are  denoted  by 

FJP  =  Jj04n*m>- 

FJPP  =  Jj<km+lRm>* 

FNP  =  N!(k  R  ),  and 
J  m  m 

FNPP  =  N!(k  ,R  )  . 

j  m+ 1  m 

Between  S80  and  SluO,  the  constants  Umn,  Vmn»  Wmn,  and  Xmn 
defined  in  Eqs.  (30)  through  (33)  are  calculated  and  are  denoted  by  UJ» 
VJ,  WJf  and  XJ  for  the  TM  case,  and  by  UP,  VP,  WP,  and  XP  for  the 
TE  case.  Equations  (28)  and  (29)  are  employed  to  calculate  A^n+j  n 
and  Bm+l,n  which  are  denoted  by  A(J)  and  B(J)  for  the  TM  case  and  by 
AP(J)  and  BP(J)  for  the  TE  case.  Between  SI  10  and  SI 50  the  scattering 
coefficients  Dn  are  calculated  by  means  of  Eq.  (45).  They  are  denoted 
by  D(I)  and  DP(I)  for  the  TM  and  TE  cases,  respectively. 

In  the  remainder  of  the  program  Eq.  (48)  is  used  to  calculate  the 
echo  width  of  the  layered  cylinder  as  a  function  of  the  angle  <j>.  The 
symbols  for  the  output  data  are  defined  by 

PHI  =  Scattering  angle  9  in  degrees, 

ECHW  =  Echo  width/ wavelength  for  TM  case,  and 

ECHWP  =  Echo  width/ wavelength  for  TE  case. 
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*#* 

#«# 


FIRST 


S25 


S26 


S50 


RUN 

SCATRAN 

COMPLEX (EXP#ES*FF#FH«D#OP#ESP )  — 

D1MENSI0N(RL( 102) #E< 102) #U(l02) #D< 102) #FJA< 102)#DP< 102 ) « 

FJBU02)  *FNA(  102)  #FNB<  1  02  )  .  A  (1  02  )  .9  ( 1 02  )  •  AP  (  1  02  )  #BP  (  1  02  )  ) 

P I  *  3# 1 4 1 5926- 

TP=2**PI- 

P2=PI/2*- 

CST*2# /P I  — 

READ  INPUT# 7* (FN#DPH)- 

n*fn- 

READ  INPUT #7. < ( RL (  I  )  •  1  =  1  •  1  •  I #LE#N) )- 
READ  INPUT# 7# ( ( E ( I ) # I = 1 . 1 • I *LE.N) )- 
READ  INPUT #7# ( <U( I ) . I = 1 # 1 • I #LE#N) >- 
DO  THROUGH ( S25) #1=1 #1#I#LE#N- 
WRITE  OUTPUT #2 ♦ (RL( I)#E(I).U(I))- 
WRITE  OUTPUT ♦  2- 
NN=N+ 1  — 

Lc4#+TP*RL<N)*SQRT# (E(N)*U(N) )- 
PROV I DED  J L#G# 1 00 ) «L=1 00- 
LL=L- 1 - 
E ( NN ) = I • — 

U(NN)=1 •  - 
I JK«L- 

DO  THROUGH (S26)«I: 0*1# I.LE.L- 
A ( I ) = 1 •- 
8 ( I ) = • 0- 
AP ( I ) = 1 • — 

BP ( I )=.0- 

DO  THROUGH (SI  10) *M»1 • 1 #M#LE*N- 
mm=m+i- 

L=I JK- 

W»TP*RL(M)*SORT. (U(M)*E(M) >- 

CALL  SUBROUT  I NECFJA«FNA#L* IND)*BFX90# (Hf)- 

WsTP*RL(M)*SORT« <U(MM)*F<MM> )- 

CALL  SUBROUT  I NECFJB.FNB.L# I ND ) =BFX90 •  ( W  )  — 

LL=L+ 1 - 

DO  THROUGH (S50 ) « IK=LL# I#IK*LE#JJK- 

FJA ( IK)=.0- 

FNA  (  IK ) s #0- 

FJB( IK)=.0- 

FNB  { IK ) * • 0— 

CNST  =  P2*U ( MM ) *TP*PL ( M ) - 
CNSTP*P2*E ( MM ) *TP*RL (M )  - 
YM*SOPT# (E(M)/U(M) )*CNST- 
YMP*SQRT • <U<M)/E(M) )*CNSTP- 
YMM*SORT# (E(MM)/U(MM ) )*CNST- 
YMMP=50RT. (U(MM )/E(MM ) )#CNSTP~ 

DO  THROUGH ( Si  00 )«J=0« 1 «J#LE  #LL - 

JJeJ+1- 

FJP*-FJA( JJ)- 

FJPP=-FJB(JJ)- 

FNP*-FNA ( JJ)- 

FNPPs-FNB ( JJ)- 

PROV I DED ( J#E • 0 ) ♦ TRANSFER  T0(S80)- 

JJJr.J-J.. 

FJP*#5*(FJP+FJA( JJJ) )- 


Fig.  6.  Digital- computer  program  for  scattering  patterns  of  a 
circular  dielectric  cylinder  having  several  layers. 
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sao 


si  oo 

si  10 


S  1 50 


S  1 80 


S200 

### 


Fig 


F JPP*.5*  (  FJPP+F JB  (  JJJ  )  >- 
FNP*.5* (FNP+FNA { JJJ  > )- 

FNPP*.5»(FNPP+FNB( JJJ ) )- 
FJZ*FJA< J>- 
F  J 1 =r JB ( J  )  - 
FNZ*FNA( J)- 
FN1 *FNB ( J  )  — 

UJ*VMM#FJZ*FNPP-YM*FJP*FN1 - 
VJ*YM#F JP»F J1 -YMM#FJZ*FJPP- 
WJ*YMM#FNZ*FNPP“YM#FNP#FN1 — 

XJ=YM*FNP*F J1 -YMM*FNZ*F JPP- 
A J= A ( J ) - 
B J=B ( J  )- 

A  ( J  )=UJ*AJ4-WJ*BJ- 
B ( J ) *V J*A J+X J*B J- 
UP*YMMP#F JZ*FNPP-YMP»F JP*FN1 - 
VP=YMP*FJP#FJ1-YMMP#FJZ*FJPP- 
WP»YMMP#FNZ*FNPP-YMP*FNP#FN1 - 
XP=YMP#FNP*FJ1 -YMMP*FNZ*FJPP- 
A J=AP{ J) - 
RJ=BP ( J ) - 

AP ( J ) *UP* A J+WP*BJ- 
BP( J ) *VP*A J+XP*BJ- 
CONTINUF- 
CONT INUE- 

DO  THROUGH (Si  50). I *0.1.1 .LE.LL- 
F 1  =  I  — 

D ( I ) =-B ( I )/(B(I)+(.I.l.)*A(I))~ 

0P( I )=-BP( I )/(BP( I  )+(• I .1 • )*AP( I ) )- 
DMAG=.ABS.D( 1 >- 
DM AGP* • ABS.DP ( I >- 

write  output . i «  cfi .d< j >  .dpci >.dmag.dmagp>- 

WRITE  OUTPUT *2- 
M=180./0PH- 

DO  THROUGH (S200 ) • I =0* 1 • I .LF.M- 
ES*  .0— 

ESP* • 0- 
F!*I- 

PH I *F I *DPH— 

PR*  *01 745320*  PH I  — 

EN* 1 •  - 

DO  THROUGH (SI  80 )«J=0«1 .J.LE.LL- 
FF J= J— 

COS*COS.  (FFJ#PR)- 
ES*ES+EN*D ( J ) *COS- 
ESP=ESP+EN#DP  ( J  )»COS- 
EN=2 • — 

EMAGP=*ABS.ESP- 

ECHWP*CST#EMAGP»EMAGP- 

EMAG*.ABS.ES- 

ECHW=CST*EMAG»EMAG- 

WRITE  OUTPUT. 2« (PHI . ECHW.ECHWP )- 

END  PROGRAM (FIRST )- 

DATA 


5.0 

•  1 

C  OJ 

.  . 

in 

.3 

•  4 

.5 

6.0 

5.0 

4.0 

3.0 

2.0 

1  .0 

1.0 

1.0 

1  .0 

1.0 

6.  Digital- c ornpute r  program  for  scattering  patterns  of 
circular  dielectric  cylinder  having  several  layers. 
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